home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / ch_3.1 / xconv / xconv.c < prev    next >
C/C++ Source or Header  |  1999-09-11  |  7KB  |  300 lines
  1. /* 
  2.  * xconv.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * XCONV(olve) 
  11.  *
  12.  * test convolve.c module
  13.  * 
  14.  */
  15.  
  16. #include "xconv.h"
  17.  
  18. #define BUFSZ 1024
  19. #define    MAX_FILTER_LENGTH    65
  20. #undef FILTER_PRINT             /* define this to print filter */
  21.  
  22. /*
  23.  * globals
  24.  */
  25. extern short tiffInput;         /* flag=0 if no ImageIn to set tags; else =1 */
  26. extern char *optarg;
  27. extern int optind, opterr;
  28.  
  29. /*
  30.  * usage of routine
  31.  */
  32. void
  33. usage (char *progname)
  34. {
  35.   progname = last_bs (progname);
  36.   printf ("USAGE: %s inimg outimg [-f filter_file] [-g len] [-L]\n", progname);
  37.   printf ("\n%s convolves an image with a kernel;\n", progname);
  38.   printf ("kernels may be read from a filter file; Gaussian\n");
  39.   printf ("kernels of any length can be generated automatically.\n\n");
  40.   printf ("ARGUMENTS:\n");
  41.   printf ("    inimg: intput image filename (TIF)\n");
  42.   printf ("   outimg: output image filename (TIF)\n\n");
  43.   printf ("OPTIONS:\n");
  44.   printf ("  -f filter_file: filename for filter to use (ASCII)\n");
  45.   printf ("          -g len: use 2-D Gaussian filter of length len\n");
  46.   printf ("                  NOTE: either -f or -g must be specified\n");
  47.   printf ("              -L: print Software License for this module\n");
  48.   exit (1);
  49. }
  50.  
  51.  
  52. #undef FILTER_TEST
  53.  
  54. Matrix *
  55. genGaussFilter (char *filterLenStr)
  56. {
  57.   Matrix *retMatrix;
  58.   float *ip, **ipp;
  59.   int ix, iy;
  60.   int i;
  61.   int nf2;
  62.   float two_sigma_sq;
  63.   float x;
  64.   float sigma;
  65.   float *filter1D;
  66.   float norm;
  67.   int filterLen;
  68.  
  69.   if (!sscanf (filterLenStr, "%d", &filterLen)) {
  70.     printf ("Can't scan gaussian filter length\n");
  71.     exit (1);
  72.   }
  73.   if (filterLen > MAX_FILTER_LENGTH || filterLen < 0) {
  74.     printf ("\n...filter length too long!! forget it!!\n");
  75.     exit (1);
  76.   }
  77.  
  78.   filterLen = (filterLen % 2 > 0 ? filterLen : (filterLen + 1));
  79.   sigma = (float) filterLen / (float) (3.0 * 2.0);
  80.   nf2 = filterLen / 2;
  81.   norm = (float) 0.0;
  82.  
  83.   if ((retMatrix = (struct Matrix *) calloc (1, sizeof (struct Matrix))) == NULL) {
  84.     printf ("\n...mem allocation for retMatrix failed\n");
  85.     exit (1);
  86.   }
  87.   if ((ipp = (float **) calloc (filterLen, sizeof (float *))) == NULL) {
  88.     printf ("\n...mem allocation for ipp failed\n");
  89.     exit (1);
  90.   }
  91.   if ((ip = (float *) calloc (filterLen * filterLen, sizeof (float))) == NULL) {
  92.     printf ("\n...mem allocation for ip failed\n");
  93.     exit (1);
  94.   }
  95.   if ((filter1D = (float *) calloc (filterLen, sizeof (float))) == NULL) {
  96.     printf ("\n...mem allocation for filter1D failed\n");
  97.     exit (1);
  98.   }
  99.   for (iy = 0; iy < filterLen; iy++)
  100.     ipp[iy] = ip + (iy * filterLen);
  101.   retMatrix->matrix = ipp;
  102.   /*
  103.    * Generate a 1-D Gaussian
  104.    */
  105.   two_sigma_sq = (float) 2.0 *sigma * sigma;
  106.   for (i = 0; i < filterLen; i++) {
  107.     x = (float) (i - filterLen / 2);
  108.     *(filter1D + i) = (float) (255.0 * exp (-(x * x) / two_sigma_sq));
  109.     norm += *(filter1D + i);
  110.   }
  111.   norm = (float) fabs ((double) norm);
  112.  
  113.   /*
  114.    * Now create the 2-D kernel
  115.    */
  116.   for (iy = 0; iy < filterLen; iy++)
  117.     for (ix = 0; ix < filterLen; ix++)
  118.       ipp[iy][ix] = (*(filter1D + iy) * *(filter1D + ix));
  119. #ifdef FILTER_PRINT
  120.   printf ("...sigma = %f\n", sigma);
  121.   printf ("...filter length = %d\n", filterLen);
  122.   printf ("...kernel coefficients:\n");
  123.   for (iy = 0; iy < filterLen; iy++) {
  124.     for (ix = 0; ix < filterLen; ix++)
  125.       printf ("%5.5f  ", ipp[iy][ix]);
  126.     printf ("\n");
  127.   }
  128. #endif
  129.   retMatrix->nRows = filterLen;
  130.   retMatrix->nCols = filterLen;
  131.   return (retMatrix);
  132. }
  133.  
  134. Matrix *
  135. readFilter (char *filename)
  136. {
  137.   FILE *stream;
  138.   int nx, ny;
  139.   Matrix *retMatrix;
  140.   float *ip, **ipp;
  141.   int ix, iy;
  142.   char *sptr, inBuf[BUFSZ];
  143.   char seps[] = " ,\t\n";
  144.   char *token;
  145.  
  146.   if ((stream = fopen (filename, "r")) != NULL) {
  147.     for (;;) {
  148.       if (fgets (inBuf, BUFSZ, stream) == (char *) NULL) {
  149.         printf ("Premature EOF encountered while reading filter file %s\n", filename);
  150.         exit (1);
  151.       }
  152.       /*
  153.        * Throw away comments
  154.        */
  155.       if (inBuf[0] == '#')
  156.         continue;
  157.       else {
  158.         if (sscanf (inBuf, "%d %d\n", &nx, &ny) != 2) {
  159.           printf ("error getting kernel size\n");
  160.           exit (1);
  161.         }
  162.         if ((retMatrix = (struct Matrix *) calloc (1, sizeof (struct Matrix))) == NULL) {
  163.           printf ("\n...mem allocation for retMatrix failed\n");
  164.           exit (1);
  165.         }
  166.         if ((ipp = (float **) calloc (ny, sizeof (float *))) == NULL) {
  167.           printf ("\n...mem allocation for ipp failed\n");
  168.           exit (1);
  169.         }
  170.         if ((ip = (float *) calloc (ny * nx, sizeof (float))) == NULL) {
  171.           printf ("\n...mem allocation for ip failed\n");
  172.           exit (1);
  173.         }
  174.         for (iy = 0; iy < ny; iy++)
  175.           ipp[iy] = ip + (iy * nx);
  176.         retMatrix->matrix = ipp;
  177.         break;
  178.       }
  179.     }
  180. /*
  181.  * Now read the kernel in
  182.  */
  183.     for (iy = 0; iy < ny; iy++) {
  184.       if (fgets (inBuf, BUFSZ, stream) == (char *) NULL) {
  185.         printf ("Premature EOF encountered while reading filter file %s\n", filename);
  186.         exit (1);
  187.       }
  188.       /*
  189.        * Throw away comments
  190.        */
  191.       if (inBuf[0] == '#')
  192.         continue;
  193.       sptr = inBuf;
  194.       for (ix = 0; ix < nx; ix++) {
  195.         if ((token = strtok (sptr, seps)) == (char *) NULL) {
  196.           printf ("Can't read filter values on line %d of file %s\n", iy + 2, filename);
  197.           exit (1);
  198.         }
  199.         if (!sscanf (token, "%f", &(ipp[iy][ix]))) {
  200.           printf ("Can't scan filter values on line %d of file %s\n", iy + 2, filename);
  201.           exit (1);
  202.         }
  203.         sptr = NULL;
  204.       }
  205.     }
  206.     fclose (stream);
  207.     retMatrix->nRows = ny;
  208.     retMatrix->nCols = nx;
  209.     for (iy = 0; iy < ny; iy++) {
  210.       for (ix = 0; ix < nx; ix++)
  211.         printf ("%5.5f  ", ipp[iy][ix]);
  212.       printf ("\n");
  213.     }
  214.     return (retMatrix);
  215.   }
  216.   else {
  217.     printf ("Can not open kernel file: %s\n", filename);
  218.     exit (1);
  219.   }
  220. }
  221.  
  222. int
  223. main (int argc, char *argv[])
  224. {
  225.  
  226.   Image *imgIn, *imgOut;
  227.   Matrix *kernelp;
  228.   int kernel_gauss = 0;
  229.   int kernel_file = 0;
  230.   int i_arg;
  231.  
  232. /* 
  233.  * cmd line options:
  234.  */
  235.   static char *optstring = "f:g:L";
  236.  
  237. /*
  238.  * parse command line
  239.  */
  240.   optind = 3;
  241.   opterr = ON;                  /* give error messages */
  242.  
  243.   if (argc < 4)
  244.     usage (argv[0]);
  245.  
  246.   while ((i_arg = getopt (argc, argv, optstring)) != EOF) {
  247.     switch (i_arg) {
  248.     case 'f':
  249.       /*
  250.        * Read in filter from file
  251.        */
  252.       kernelp = readFilter (optarg);
  253.       kernel_file = 1;
  254.       break;
  255.     case 'g':
  256.       /*
  257.        * Generate 2-D Gaussian
  258.        */
  259.       kernelp = genGaussFilter (optarg);
  260.       kernel_gauss = 1;
  261.       break;
  262.     case 'L':
  263.       print_sos_lic ();
  264.       exit (0);
  265.     default:
  266.       usage (argv[0]);
  267.       break;
  268.     }
  269.   }
  270.   if (kernel_gauss && kernel_file) {
  271.     printf ("Please specify either -f OR -g option, not both!\n");
  272.     exit (1);
  273.   }
  274. /*
  275.  * read input image
  276.  */
  277.   imgIn = ImageIn (argv[1]);
  278.   if (imgIn->bps == 8 && imgIn->spp == 3) {
  279.     printf ("Got RGB image!!!\nInput image must be Grayscale or B&W!!\n");
  280.     exit (1);
  281.   }
  282.   /* reset tiffInput so that we write a grayscale file (i.e tags are not copied) */
  283.   tiffInput = 0;
  284. /*
  285.  * Allocate memory for output image
  286.  */
  287.   imgOut = ImageAlloc (imgIn->height, imgIn->width, imgIn->bps);
  288.  
  289.   printf ("Doing image convolution...\n");
  290. /*
  291.  * Do the convolution
  292.  */
  293.   convolve (imgIn, imgOut, kernelp);
  294. /*
  295.  * Write the output image
  296.  */
  297.   ImageOut (argv[2], imgOut);
  298.   return (1);
  299. }
  300.